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Using a recently presented numerical code for calculating the Lorenz-gauge gravitational self-force 
(GSF), we compute the 0(m) conservative correction to the precession rate of the small-eccentricity 
orbits of a particle of mass m moving around a Schwarzschild black hole of mass M S> m. Specifically, 
we study the gauge-invariant function p(x), where p is defined as the 0(m) part of the dimensionless 
ratio (Cl r /ti v ) 2 between the squares of the radial and azimuthal frequencies of the orbit, and where 
x = [Gc~ 3 (M + m)f2 v ] 2//3 is a gauge-invariant measure of the dimensionless gravitational potential 
(mass over radius) associated with the mean circular orbit. Our GSF computation of the function 
p{x) in the interval < x < 1/6 determines, for the first time, the strong-field behavior of a 
combination of two of the basic functions entering the Effective One Body (EOB) description of 
j^Q 1 the conservative dynamics of binary systems. We show that our results agree well in the weak-field 

r-{ ' regime (small x) with the 3rd post-Newtonian (PN) expansion of the EOB results, and that this 

agreement is improved when taking into account the analytic values of some of the logarithmic- 
running terms occurring at higher PN orders. Furthermore, we demonstrate that GSF data give 
' access to higher-order PN terms of p{x) and can be used to set useful new constraints on the 

values of yet-undetermined EOB parameters. Most significantly, we observe that an excellent global 
representation of p(x) can be obtained using a simple 'two-point' Pade approximant which combines 
3PN knowledge at x — with GSF information at a single strong- field point (say, x = 1/6). 
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I. INTRODUCTION 



The problem of calculating the gravitational self-force (GSF) acting on a pointlike particle of mass m, in bound 
orbit around a Schwarzschild black hole (of a much larger mass M), is now well understood at leading order in the 
£f~) | small mass ratio q = m/M (see 0, Q for reviews of recent developments). Explicit numerical computations of the 
Q\ . GSF, specialized to circular orbits, were carried out by at least three different groups using a variety of methods 

These calculations all rely on mode-sum regularization [(| 0] but are otherwise based on very different computational 
strategies which are formulated in different gauges. Two of us recently reported a numerical code for computing the 
f^i 1 GSF along generic (bound) geodesies with arbitrary eccentricities ||. 

In principle, knowledge of the GSF (along with the metric perturbation associated with the particle) gives a 
complete information about the O(q) post-geodesic dynamics of the binary system. However, the GSF itself is a 
gauge- dependent notion [9(, and it is generally not a straightforward task to translate the raw GSF data (i.e., the 
value of the various components of the GSF along the orbit, computed in a given gauge) into physically meaningful 
gauge-invariant statements about the orbital dynamics. It is particularly non-trivial to obtain such a gauge-invariant 
description for the conservative part of the dynamics, which is also (in principle) accessible from the GSF. For circular 
orbits, Detweiler Q proposed to measure the conservative effect of the GSF using the gauge-invariant relation between 
the "red-shift" parameter u (the Schwarzschild-i contravariant component of the four-velocity along the circular 
geodesic) and the azimuthal frequency Ct v . Using numerical data for the GSF in the Regge- Wheeler gauge, Detweiler 
was able to compute the relation through 0(q). Numerical computations of the GSF in the Lorenz-gauge and 

in the radiation gauge were later shown to reproduce the same relation u*(f2 v ) fiol. ITTI]. 

The recent development of a GSF code for eccentric orbits (formulated in the Lorenz gauge) enabled the computation 
of additional gauge invariant quantities associated with the conservative effect of the GSF. By considering slightly- 
eccentric orbits, two of the authors have computed the O(q) conservative shift in the orbital frequency of the innermost 
stable circular orbit (ISCO) H, E3- This quantity (unlike it*) has a clear physical interpretation, and, at least in 
principle, a measurable effect on the signature of gravitational waves from the binary system. Most recently, the GSF 
correction to the precession rate of orbits with arbitrary eccentricities has also been calculated [l3j . 

Analyzing the gauge-invariant effects of the GSF allows one not only to compare between GSF calculations carried 
out in different gauges, but also to make an important connection with the well-established post-Newtonian (PN) 
theory. Indeed, in Ref. Q Detweiler used the relation w*(f2 v ) to demonstrate, for the first time, that conservative 
GSF results (for circular orbits) were consistent with the analytic predictions of PN theory. This comparison was 
further improved, and pushed to higher PN orders, in Refs. [l4j and [TBj . 
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Such "cross-cultural" comparisons (to use the language of Ref. [14| ) are motivated in several ways, (i) They provide 
highly non-trivial opportunities to confirm some of the basic aspects of both GSF and PN formalisms (such as the 
regularization procedures underpinning them), (ii) Comparison of concrete results from GSF and PN computations 
serves to test both sets of results; such two-way validation tests are crucial given the complexity of both GSF numerical 
codes and PN analytic computations, (iii) Results from GSF calculations can inform the development of accurate 
Analytical Relativity (AR) models by providing high-accuracy calibration data for yet-unknown high-order PN terms. 
Combined also with Numerical Relativity (NR) results for comparable-mass binaries, one eventually hopes to develop 
a reliable and accurate AR model across the full range of binary mass ratios (e.g., for gravitational- wave detection 
applications) . Within such a program, GSF calculations provide a crucial "data point" at the extreme edge (extreme 
mass ratio, small separation) of the essential parameter space of the binary problem, (iv) Since PN theory is formulated 
for arbitrary mass ratios, it could be used to test second- (and higher-)order GSF computations when these become 
available in the future. 

A fundamental challenge in any GSF-PN comparison arises from the fact that the two methods are designed 
to explore distinct regimes of the binary problem: GSF methods work best for strong-field orbits (the necessary 
computational resources tend to increase fast with the orbital radius Q), while the standard PN series is essentially 
an asymptotic expansion about the limit of infinite binary separation — although it performs surprisingly well at rather 
small separations, it ultimately fails to converge at sufficiently small radii. The upshot is that, in order to facilitate 
a useful GSF-PN synergy, one must (generally speaking) stretch GSF techniques beyond their natural domain of 
applicability (i.e., obtain large-radius data) while at the same time attempt to compute new higher-order terms in 
the PN series, at least through O(q). That this can be achieved in practice was indeed demonstrated by Blanchet et 
al. in Refs. 0,111. 

Although PN theory is the basis of most of our current analytical knowledge of the motion and radiation of 
gravitationally interacting comparable-mass binary systems, it is unable, by itself, to provide an analytical description 
of the entire evolution of such systems. Indeed, the PN description is only accurate during the early inspiral and it 
breaks down during the late inspiral, well before the merger. The Effective One Body (EOB) formalism (la - [l9| was 
proposed as a flexible AR framework for describing the motion and radiation of coalescing binaries over the entire 
merger process, from the early inspiral, right across the innermost stable orbit, through the 'plunge' and to the merger 
and final ringdown. 

At the heart of the EOB formulation is an "effective Hamiltonian" which depends on 3 initially-unspecified functions: 
A{u\v), D{u\v) and Q{u,p r ;v) (see Ref. [2(| for a review). Here u = G(M + m) / (c 2 t-eob) is the dimensionless 
gravitational potential (with teob being the EOB measure of the binary separation), v = mM/(m + M) 2 is the 
symmetric mass ratio, and p r is the EOB (relative) radial momentum. PN theory gives access to the expansions of 
these functions in powers of the variable u. These "Taylor-expanded" functions were found to behave badly (at the 
3PN level) in the strong- field region u ~ 1/6 [lH. It was then suggested to replace them with suitably resummed 
expressions (Pade approximants). The natural analytical flexibility of the EOB formalism leads one to consider 
resummed functions, especially A(u; v\ which include yet uncalculated coefficients, as,a6, ... parametrizing PN terms 
beyond the currently known 3PN level. These a priori unknown coefficients can then be "calibrated" by comparing 
EOB predictions to numerical data from NR simulations (2l| - |24| . Vigorous ongoing activity in this area is convincingly 
demonstrating the utility of the EOB framework to describe accurately all phases of the merger process of comparable- 
mass binaries, from the early inspiral to the final ringdown. Moreover, EOB theory has been shown to perform very 
well also in the extreme mass ratio regime [25l - l27l |. 

The recent GSF results provide a new source of calibration data for the EOB model. The GSF data have the 
potential to be particularly useful for this purpose, for a number of reasons, (i) Unlike most available NR data, GSF 
data is "clean" and very accurate — the numerical component in GSF calculations only involves linear differential 
equations, (ii) In GSF calculations, again unlike in NR, it is straightforward to disentangle the conservative aspects 
of the dynamics from the dissipative ones (EOB treats these two aspects separately), (iii) GSF data "fill the gap" 
in the strong-field/extreme- mass-ratio corner of the parameter space, which (as of yet) is inaccessible to either PN 
or NR. (iv) GSF calculations already return accurate data for orbits with large eccentricities; this data could give a 
good handle on the two EOB functions {D and Q) which describe the radial component of the motion and have not 
been studied in detail so far. On the other hand, the GSF data have the shortcoming of only giving access to the 
0{v) terms in the ^-expansions of the EOB functions A{u] v), D(u; v) and Q(u,p r ; v). 

The prospect for a useful GSF-EOB synergy was recentlyhighlighted in a paper by one of us [Hj]. This work 
explored how the single GSF data point obtained in Ref. [12[ [namely, the 0(q) conservative shift in the ISCO 
frequency] constrains the shape around u = 1/6 of the 0(v) term in the crucial EOB radial potential A(u;v), and 
thereby can improve the determination of two higher-order parameters 05 and ci6, which previous EOB-NR analysis 
had found to be strongly degenerate. Ref. [28| also identified several additional gauge-invariant quantities that might 
be computable using current GSF technology, and could provide useful additional calibration data for EOB. 

In this work we consider one of the gauge-invariant quantities proposed in Ref. 28] , namely the conservative GSF 
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correction to the periastron precession for slightly eccentric orbits. More precisely, we consider the related quantity 
p{x), denned as the 0{q) conservative part of the squared ratio of the radial and azimuthal frequencies [see Eq. (Til))) 
below], with x being the dimensionless gravitational potential defined from the invariant azimuthal frequency. The 
function p{x) is a truly dynamical gauge-invariant characteristic of the conservative dynamics, in that, as shown in 
(28|, its knowledge gives a direct handle on the strong-field behavior of [the 0(v) part of] a combination of the EOB 
functions A(u] v) and D(u;u) . The GSF computation of the function p(x) in the interval < x < 1/6, which we 
present here, is therefore the first determination of the strong-field behavior of a combination of functions which are 



of direct significance for describing the conservative dynamics of binary systems 37[ • As such it can contribute to the 
development of accurate AR models of the gravitational wave signals emitted b y co alescing binaries. 

The weak-field behavior of p(x), on the other hand, has already been analyzed 28, 29] through 3PN order [/9 PN (.t) = 
P2X 2 + P3X 3 + 0(x 4 lnx)]. In addition, the logarithmic contributions present at the 4PN and 5PN levels (p^x 4 lnx + 
• • • + p)? 6 x 5 lnx + • • •) have been recently determined analytically [3(|. [The presence of logarithmic contributions 
in the (near-zone) PN expansions, starting at the 4PN level, follows from old work [3lT-33]. The importance and 
detectability of these logarithmic terms in the comparison between GSF data and PN expansions has been recently 
discussed in Refs. [TU, l28j].] Here we shall use the numerical methods of Refs. [1, HH to obtain corresponding GSF 
data for p(x), both in the strong-field region (x ~ 1/6) and in the weak field (x <C 1), and explore what can be learned 
by comparing the GSF results with the EOB/PN expressions. 

The paper is structured as follows. In the next section we derive an expression for the function p{x) in terms of 
GSF quantities, and review the EOB expression for p{x). In Sec. Ill we describe the numerical method applied to 
obtain the necessary GSF data, and present our numerical results for p(x). In Sec. IV we test the numerical data 
against the currently available EOB/PN results, and in Sec. V we examine to what extent the GSF data can be used 
to determine yet-unknown higher-order PN terms in the EOB model. Section VI explores the utility of simple Padc 
models (i.e., rational- function fits), based on a minimal amount of EOB and GSF information, to provide accurate 
global fits for p{x). Section VII contains a summary and a discussion of future directions for EOB-GSF synergy. 

Throughout this paper we use standard geometrized units (with G = c = 1), metric signature — I — | — h, and 
Schwarzschild coordinates (i, ?*, 6, ip). 



II. SMALL-ECCENTRICITY PRECESSION EFFECT IN THE GSF AND EOB FORMULATIONS 

We consider a gravitationally-bound binary comprising a particle of mass m and a Schwarzschild black hole of a 
much greater mass, M >m. We denote the small mass ratio by q = 77i/M(^ 1), and throughout our analysis work 
through 0(q) only. Readers should be wary of the conflict between the standard PN/EOB notation and the common 
GSF one (see Table [I]); our notation here represents a compromise between the two sets of notations. Note that 
the symmetric mass ratio [38[ is v = q/(l + q) 2 = q + 0(q 2 ), so that an O(q) quantity is also 0{y). One should, 
however, beware of the important O(q) difference between the large mass M (which is often used in GSF works to 
adimcnsionalizc frequencies, and denoted there M) and the total mass M + m = M(l + q) (which is used in PN/EOB 
works to adimensionalize frequencies, and denoted there M as well). 





This paper 


GSF literature (e.g., [8, 10]) 


PN/EOB literature (e.g., [28]) 


particle mass 


m 


P 


mi 


black hole mass 


M 


M 




total mass 


M + m 




M = mi + mi 


"small" mass ratio 


q = m/M 






symmetric mass ratio 


v = mM/(m + M) 2 




v = m 1 m2/(m,i + m^) 2 


effective mass 


a = mM/(m + M) 




fj, = mim 2 /(mi + 7-712) 



TABLE I: Our notation for various mass quantities, compared with the common GSF notation and with the standard PN or 
EOB notation. Readers should be wary of these notation differences, which can bring confusion. 
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A. GSF treatment 

1. GSF-corrected circular orbits 

In GSF theory the binary dynamics is described perturbatively: At the limit to — > the orbit is a geodesic of the 
"background" Schwarzschild geometry of mass M, and we consider the O(q) perturbation to this geodesic due to the 
effect of the GSF. Here we focus on the conservative effect alone, and ignore dissipation. We start by considering 
circular orbits, parameterized by their (Schwarzschild) radius r = r = const. (Here and in the following, ro denotes 
the radius of the perturbed, i.e., GSF-corrected circular orbit. Such orbits exist because we are considering here 
only the conservative part of the GSF.) Without loss of generality we let the orbit lie in the equatorial plane. The 
Schwarzschild components of the particle's (GSF-corrected) four- velocity u a = dx a jdr (r denoting the background 
proper time along the perturbed orbit) can be written as 

«« = ^{1,0,0,11,}, (i) 

where E = ~-gt a u a is the specific (background-defined) 'energy' of the particle, and Q v = dip/dt is the azimuthal 
coordinate-time (tp-) frequency. The perturbed values of the latter two quantities (squared, for convenience) are related, 
through 0(g), to the radius r of the orbit by [3| 



e2= (TTQ-2M) 



r (r -3M) 



r, 



1 y p r 



'o 



m(r - 2M)' 
rg(r -3M) 



circ 



mM(r - 2M) 



F r - 



(2) 



(3) 



where F c r irc ((x to 2 ) is the radial component of the GSF acting on the particle (which, in the case of circular motion, 
is entirely conservative). We note that, in practice, it is sufficient for us to evaluate the GSF along the background 
geodesic rather than along the perturbed orbit, as the resulting difference in the value of dynamical quantities like E 
or Q v is only 0(g 2 ) and can be neglected here. 

We recall that in GSF theory the radius ro is gauge dependent, just like the GSF itself |9]. In this work we will 
always consider the GSF in the Lorenz gauge, which is the choice of gauge made in Refs. [H, HI]. The GSF component 
F£ irc , and all other GSF quantities introduced below, should be understood to be given specifically in the Lorenz 
gauge. We also recall that in GSF theory the value of the particle's specific energy E, given in Eq. ([2]), is dependent 
upon the gauge. However, the frequency Q v in Eq. (|3|) is gauge invariant. More precisely, Q v is invariant under 
the restricted class of 0(g) gauge transformations whose displacement vectors respect the helical symmetry of the 
perturbed spacetime (see, e.g., |10f for a detailed discussion of this po int) . 

Importantly, however — as discussed previously in the literature [13, HH, HH — the gauge transformation relating the 
Lorenz-gauge metric perturbation to the ones used in PN and EOB studies, does not fall in the above category, and 
therefore does not leave f2 v invariant. Nonetheless, it is straightforward to account for this gauge difference using 
what amounts to a simple O(g) "rescaling" of the time t (Toj . 

£—*•£= (1 + qa)t, (4) 

with 

a= MM7-0-3M)]- 1 / 2 . (5) 

This leads to a "rescaled-f frequency, given through 0(q) by 

&<p = (1 - qa)Q v . (6) 

The frequency Cl v refers to an asymptotically-flat coordinate system (as the one employed in PN and EOB studies) 
and it therefore provides a useful reference point for comparison between GSF and PN/EOB calculations. To make 
further contact with standard PN/EOB notions, we also introduce the dimensionless gravitational potential x ("total 
mass over radius"), defined through 

x= [(M + m)Sl v ] 2 / 3 . (7) 
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The quantity x is the primary gauge-invariant characteristic of the conservative dynamics of GSF-corrected circular 
orbits. Using Eq. in conjunction with Eqs. ([3]), ([5]) and ©, one finds that the relation between r and x is given 
through 0(q) by 



ro = — 

x 



1 + §, 



1 



(l-Sa;) 1 ^; 3x 2 (l 



1 



3^ pr 



(8) 



2. GSF-corrected slightly-eccentric orbits 

A second gauge-invariant quantity associated with the conservative dynamics can be constructed by considering 
a small-eccentricity perturbation of the circular orbit. Through linear order in the eccentricity e, the radius of the 
slightly-eccentric orbit can be written in the form r(r) = ro(l — ecoso; r T), where w r is the radial frequency defined 
with respect to the proper-time r along the orbit (and where, without loss of generality, we have set the orbital 
phase so that r = corresponds to a periastron passage). Henceforth, ro will denote the average radius of such a 
slightly-eccentric orbit, and we shall associate the gauge invariants Cl v and x introduced above to the 'mean' circular 
orbit corresponding to the eccentric orbit in question. 

As shown in Ref. [iiq . the various non-zero components of the conservative GSF along a slightly-eccentric orbit 
have the form [through 0(e)] 

F r = F^ rc + eFlcosuj r T, (9) 
Ft — euj r F t \ sina; r r, (10) 
F v — euj r F v ishiuj r T. (11) 

An expression for uj r at the limit e — > and through 0(q) was obtained in Ref. [12j. It reads 

2 M(r -6M) 3(r -4M) 1 2[M(r - 3M)] 1 / 2 



r rg(r -3M) mr (r - 3M) clrc mr„ 

The first term on the right-hand side of this expression would give the value of the radial frequency in the test-particle 
limit q — > 0; the subsequent terms [each of 0(q)] describe conservative GSF corrections. For our current analysis we 
also introduce the radial frequency fi r defined with respect to coordinate time t, along with its "rescaled-i" version 
Cl r . The latter is related to uj r (in the limit e — > 0) through 



- dr dt 

' l r = ~ Ul r 

dt dt 



1 - 2M/r 



E(l + qa) 

Substituting from Eqs. @ and p^|) we obtain, at the limit e — > and through O(q) 



U! r . (13) 



Sro^lOM Fr 2[M(r -3M)] 1 / 2 F 

r -2M + r 3 



(14) 



Since both frequencies fl v and f2 r are gauge invariant (in the aforementioned sense), the functional relation Ct r (Cl v ), 
or equivalently Cl r (x), is a genuinely gauge-invariant characteristic of the conservative dynamics. Following Ref. [28| . 
we introduce here the more convenient (dimensionless) quantity 

w=(n r /n v y, (15) 

and consider the equivalent gauge- invariant relation W(x). With Eqs. ([3]), dU), (JU) and (fl"4|). this relation takes the 
form 

W(x) = l-6x + qp(x) + 0{q 2 ), (16) 

where the 0(q) part is given by 

p{x) = f r0 (x)K iTC + f rl {x)F{ + f vl (x)F vl + f (a) (x). (17) 
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Here F c r irc = q F£ llc , F{ = q ^F^ and F vl = m F v i, and the various ^-dependent coefficients read 

fro{x) x*(l-2x) ' (18) 

frM = (19) 

f vl (x) = -2^(1 -3a:) 3 / 2 , (20) 
4x 2 

/(-)(*) = 4 -(T3^- ( 21 ) 

Provided with GSF data for slightly-eccentric orbits, Eq. (|17p can be used to compute the GSF-induced shift in the 
invariant function W(x). Let us note that this function has a simple physical interpretation: The quantity 

k = n v /n r - i = w~ 1/2 - i (22) 

describes the fractional periastron advance per radial period. The 0{q) GSF correction to k is given by 

5k= --{l-6x)- 3 ^ 2 qp{x). (23) 

Note that, in the test-mass limit, we have the well known ISCO behavior k ~ (1 — 6x) -1 / 2 . This singular behavior 
is avoided by working (as we do here) with the function W(x) = (1 + fc)~ 2 , which is smooth across the ISCO (and 
vanishes there) ; cf. [18|, (29| . 

It should also be noted that the quantity p(l/6) is related in a simple way to the ISCO frequency shift computed 
in Ref. [Hj]. The ISCO location is defined through W(x isco ) = 0, which, recalling Eq. (fT5|) . gives qp(l/6) = 6xi SCO — 1 
[neglecting terms of 0(q 2 )}. Using Eqs. ([7]) and (|6]) this then leads to 

p(l/6) = | [(7 _1 Ar> isco /fi isco + 1 - 1/718] , (24) 
where Af2i sco /Sli sco is the fractional O(q) shift in the (Lorenz-gauge) fl v at the ISCO. The numerical value obtained 



m 



12J for the latter was (0.4870 ± 0.0006)g, giving p(l/6) = 0.8342 ± 0.0004. As part of our current analysis we will 



obtain the more accurate value p(l/6) = 0.83413 ± 0.00004 (cf. Table HI below). 



B. EOB treatment 



An EOB treatment of slightly-eccentric orbits in Schwarzschild in the small mass-ratio case was presented in the 
recent work [28| . Here we shall merely summarize the relevant results, and we refer to reader to Ref. [28| for full 
details. 

In the EOB approach the two components of the binary system are treated "on equal footing" : Unlike in the GSF 
approach, where one speaks of a "background" and a "perturbation" , here the conservative dynamics is described 
in terms of an effective geometry which depends on the binary masses only through the symmetric combinations 
/j = Mm/(M + m) (effective mass) and v = Mm/(M + to) 2 (symmetric mass ratio; cf. Table [XJ) . The effective 
metric g^, is spherically symmetric, and is a 'deformed' version of the Schwarschild metric (of mass M + m), with 
v being the deformation parameter. When using Schwarschild-like EOB coordinates, is fully described by two 
functions of the EOB radial coordinate teob, namely A{u:v) = — g^ and D{u;v) = — (SooSrr)^ 1 - Here, the 
argument u = (M + to) /teob is the dimensionlcss EOB gravitational potential, which conveniently parametrizes the 
binary separation in the EOB formalism. As already mentioned above, the full EOB description of the conservative 
dynamics involves, besides A(u; u) and D(u; v), a third function, which depends not only on u and v but also on the 
radial (relative) momentum p,., namely Q(u,p r ; v). However, as explained in Ref. [281 ] . the conservative EOB dynamics 
of small- eccentricity orbits depends only on the two functions A{u\ v) and D(u; u) parametrizing the effective metric 

To make contact with GSF theory, one formally expands the functions A{u\ y) and D{u\v) in powers of the 
symmetric mass ratio through 0(v) (noting v = q through this order). This gives [281 ] 

A{u,v) = 1 - 2u + va(u) + 0(v 2 ) 

= l-2u + qa(u) + 0{q 2 ) 1 (25) 

D(u,v) = 1 + vd(u) + 0{v 2 ) 

= l + qd(u) + 0(q 2 ), (26) 
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where the functions a(u) and d(u) carry, in principle, all information about the conservative GSF along slightly 
eccentric orbits. [The GSF dynamics of large-eccentricity orbits would also involve a third function, say q(u,p r ).] In 
particular, it was shown in Ref. [28[ that the quantity p(x) defined in Eq. (fl6|) can be constructed from these two 
EOB functions using the formula 

p(x) = pe(x) + p a (x) + pd(x), (27) 

where 

^-^-Tnc)- (28) 

p a (x) = a(x) + xa'(x) + -x(l - 2x)a"(x), (29) 

p d (x) = (l-6x)d(x), (30) 

with a prime denoting d/dx. [Although we shall not need it here, we note that the value of the EOB gravitational 
potential u along circular orbits is related to the gauge- invariant frequency parameter x introduced above via a relation 
of the type u — x + 0(v); the 0{v) difference between u and x is given explicitly in Eq. (4.22) of [28].] Later in our 
discussion we shall occasionally refer to the 'corrected' p function obtained by subtracting from it the contribution 
Pe{x). We shall then denote it as 

p(x) = p(x) - p E {x) = p a (x) + p d (x). (31) 

In summary, a GSF computation of the function p{x) gives a direct access to a linear combination of the 0{v) EOB 
functions a(x) and d{x), and their derivatives. 



1. PN expansion 

While, as we just saw, GSF theory can give a handle on the strong-field [3!| behavior of the EOB functions, PN 
theory gives access to the weak-field behavior of the functions a(u) and d(u), i.e., to their expansions in powers of u. 
The structure of the PN expansions of these two functions reads 

a PN ( u ) = Y, a n u n , (F N (u) = J2 dnn n , (32) 

n>3 n>2 

where the (2PN and 3PN) coefficients a 3 and 04, as well as d^ and J3, are pure numbers, while higher-order coefficients 
generally involve a logarithmic dependence on u. Inserting the expansions (|32[) (with the replacement u — > x (40| ) 
into Eqs. (f27 |) - ((3"0"j) yields the weak-field (x — > 0), or PN, expansion of the function p(x), say 

P ™{x) = Y,PnX n . (33) 

ri>2 

Note that a term of order 0{x n ) in this expansion corresponds to the nth PN order [while the nPN level corresponded 
to a term of order 0(u n+1 ) in the expansion of the function a(u)]. The coefficients of the first two terms in the 
expansion p PN (x), corresponding to the 2PN and 3PN levels, are pure numbers, while the higher-order coefficients 
p n , for n = 4, 5, . . . inherit a logarithmic dependence on x from that of 05, as, . . . and d^, d§, . . .. We shall then write 
the PN expansion of p{x) in the explicit form 

p™( x ) = P2X 2 + p ^ + (pi + p'°s\nx)x i + {pi + pi° s \nx)x 5 + O(x 6+0 ), (34) 

where p\, p l £ s , p\ and p l £ s are all numerical coefficients, and where the symbol O(x 6+0 ) refers to the presence of 
logarithmic corrections (of a finite, but unspecified, order) in the 0(x 6 ) remainder. [As argued in [TBI, [28| one expects 
no higher powers of In a; to occur through 0(x 5 ), although such terms may well appear at higher orders.] 

The first two terms in the PN expansion (1341) are determined by the currently known 3PN results. More precisely, 
inserting the 3PN results [l8| 

a 3 - 2, (35) 
94 4l7r 2 

a 4 = y - -gg- 18.687903, (36) 

d 2 = 6, (37) 
d 3 = 52 (38) 
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into the first two equations obtained by inserting the expansions (|32[) into Eqs. ([27)) - (l30l) . namely 

p 2 = 2 + 3a 3 + J 2 , (39) 

P3 = -~- 2(23 + 6(14-6(22 + 4, (40) 

yields 

P2 = 14, 

397 123 9 

p 3 = 7T 2 ~ 122.627416. 41) 

H 2 16 v ; 

Recently, the logarithmic contributions to p(x) at the 4PN and 5PN levels, i.e., the analytic values of p l £ s and p l £ s , 
have been derived by one of us [30(, using effective-action techniques, with the results 

16 

p 4 og = — x 157 ~ 167.466666, 
15 

l%* = ~ -1619.428571. (42) 

By contrast, the values of the non-logarithmic coefficients p\ and p§ remain unknown, being beyond the present 
capabilities of PN theory. 

A GSF computation of p{x) can be tested against the PN expression (|34[) (with the known values of p2, p3, p l ± s 
and P5° s ), and, furthermore, can in principle be used to determine p\, p|, and possibly some other, higher-order terms 
in the expansion. For this purpose, it is useful to examine the PN behavior of the GSF expression for p(x), Eq. (|17p. 
Expanding the coefficients in Eqs. (|T8")) -(f2"l" j) in powers of x, we have 



p(x) = (~2x~ 2 + ix^ 1 + 2 + 4x + 8x 2 + . . .) F c 



circ 



+ (x~ 2 - 3a;- 1 ) F{ 

+x 1/2 (-2 + 9x - 27x 2 /A + ...) F vl 

+4x -4x 2 -6x 3 + .... (43) 

Note that, when x — > 0, Eq. (f3~4"|) says that p(x) vanishes proportionally to x 2 , while the first two terms in the 
GSF expression Eq. (|4"3")l involve coefficients that blow up proportionally to x~ 2 . This means that there must exist 
delicate cancellations between the various terms on the right-hand-side of Eq. (|43|) . It is well-known (see, e.g., (35[) 
that the leading weak-field contribution to the GSF is of Newtonian origin, and simply comes from the 'recoil' of 
the large mass with respect to the center of mass of the binary system. This leads to a leading-order GSF given by 
p(leadmg) _ +2m 2 r/r 3 , i.e., a purely radial GSF, which [upon inserting r(r) = r (l — e cos w,.r)] yields -F c r i r 1 c eadms ^ = 2x 2 
and jr|"( loadin s) = 2F^ eidlng ' > = 4x 2 . Inserting these results in the above expression for p(x), one finds that the first line 
on the right-hand side of Eq. (|4"3")l contributes — 4 + 0(x), while the second line contributes +4 + 0(x). As expected we 
have a cancellation of the leading-order terms, but this cancellation could leave a contribution to p(x) of order 0(x), 
i.e., of 1PN. To show analytically that the 1PN terms also cancel out in p(x) would require a lPN-accurate analytic 
expression for the GSF in the Lorenz gauge. Such analytical knowledge is unfortunately not yet available, especially 
for generic, eccentric orbits. (For circular orbits, the first few terms in the PN expansion of the Lorenz-gauge GSF 
have been estimated from a fit to numerical data — see Eq. (56) in Q.) That the above cancellation at 1PN does 
actually occur will be verified numerically below. 



III. NUMERICAL METHOD AND RESULTS 



In this section we describe the numerical computation of p(x) using GSF methods, and present the raw data coming 
out of this calculation. In subsequent sections we will analyze these data and explore what can be learned from a 
comparison with EOB predictions. 

We recall, noting Eq. (|T71) . that a computation of p{x) requires three pieces of input, namely the GSF coefficients 
^circ an d Fu>i- Recall also that F£ ilc is related to the radial (conservative) component of the GSF for a strictly 
circular orbit, while F[ and Fm\ are related to the small-eccentricity perturbations in the radial and azimuthal 
components of the conservative GSF. Section V.B of Ref. [8| described two independent strategies for determining 
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these coefficients. In the first strategy (dubbed method T) one considers a sequence of geodesic orbits with decreasing 
eccentricities, which approach the desired circular orbit along a suitable "track" in the 2-dimensional parameter 
space of eccentric geodesies. One then computes the GSF along each orbit in the sequence, and the necessary GSF 
coefficients are obtained through extrapolation to zero eccentricity (see [§[ for details). In the second strategy (method 
IT) the time-domain field equations themselves are expanded in the eccentricity parameter e through O(e), with the 
coefficient -F c r irc then obtained from the 0(e°) set of equations and the coefficients F[ and F v \ obtained from the 
0(e 1 ) ones. While method II is somewhat more difficult to implement, it is also significantly more computationally 
efficient, and allows one to obtain the necessary GSF coefficients with a greater accuracy. In Ref. [8[ (as part of the 
ISCO shift analysis) both methods were implemented in order to obtain the three GSF coefficients at the single radius 
r = 6M (x = 1/6). In that computation, method II was incorporated to obtain a high-accuracy result, with method 
/used to confirm that result. 

Here we shall use method II in order to obtain the GSF coefficients — and thereby p — for a large sample of x values 
between the ISCO (x = 1/6) and the weak-field radius ro = 80M (x — 0.0125). The precise implementation procedure 
follows closely that of Ref. [1[ and is based on the time-domain Lorenz-gauge GSF code described therein. We list 
below several minor details in which our current analysis deviates from that of Ref. @ . 

First, we note that the computation of the GSF coefficients at the ISCO in Ref. Q involved a certain extrapolation 
procedure even within method II. This is because the source terms in the 0(e) field equations become divergent as the 
radial frequency of the orbit approached zero [cf. Eqs. (F3)-(F12) in Appendix F of H], which makes it impractical 
to solve the e-perturbed equations at the ISCO itself. This additional computational burden is spared from us here: 
for each value x < 1/6 we need only solve the set of e-perturbed field equations once, for a particle at precisely the 
desired orbit. However, we also reproduce here (with better accuracy) the ISCO values of the GSF coefficients, and 
for this we use an extrapolation procedure similar to that employed in Ref. Q. 

On the other hand, our task here is made more challenging by the need to consider relatively large orbital radii 
(these are particularly interesting for the purpose of comparing with PN results), because the computational cost of 
our time-domain evolution tends to grow fast with increasing radius. To understand the reason, we note that in time- 
domain computations such as that of Ref. Q one does not (and usually can not) impose accurate initial conditions 
for the numerical evolution; instead, one simply let any spurious waves resulting from the imperfection of the initial 
data "dissipate away" over time, making sure that the time-evolution proceeds long enough for the magnitude of 
these transient waves to fall below a set threshold. The initial, transient stage of the evolution is then discarded, 
and one records only the physical, stationary late-time solution. The "relaxation" time of the spurious transient thus 
dictates the necessary evolution time (and we note that in our particular 1+1-dimcnsional implementation, the actual 
computational time is quadratic in the evolution time). Unfortunately, initial spurious waves from larger radii take 
longer to die off [41], and hence necessitate a longer evolution. At large orbital radii (x < 1/50) this becomes the 
main limiting factor in our numerics. 

In our numerical implementation we have adjusted the evolution time as a function of x in order to ensure that the 
error from residual initial waves is kept below (or at most comparable to) other sources of numerical error (we refer 
the reader to Ref. [8] for a detailed discussion of the various sources of error in our computation and of how these are 
monitored and controlled). Fortunately, it is the case that higher multipole modes of the metric perturbation decay 
faster than lower ones, which is easily exploitable in a mode-sum treatment like ours: It allowed us to reduce the 
evolution time for higher multipole numbers without compromising the numerical accuracy. Through experimentation, 
we arrived at the following practical scheme for determining the evolution time: (i) For 1/20 < x < 1/6 evolve the 
multipoles I = 0-3 for t = 1000M and higher multipoles for t = 500M. (ii) For 1/50 < x < 1/20 evolve the multipoles 
/ = 0-5 for t = 1500M and higher multipoles for t = 500M. (iii) For x < 1/50 evolve the multipoles I = 0-5 for 2000M 
and higher multipoles for t = 500M. 

Table [TTI summarizes our numerical results for p(x), and in Figure[T]we plot these data (dividing p by x 2 for clarity). 
Figure Q] also displays, for comparison, the various PN approximations p PN (x) derived from Eq. (IMl) . The agreement 
between the numerical data and the PN results is made more evident in Figure [2] where we plot the relative differences 
|p PN — p\ J p as functions of x. 

Inspecting Figures [T] and [H we can make the following preliminary observations, (i) The numerical GSF data seems 
to be perfectly consistent with the analytic PN expressions, (ii) At sufficiently small x (large ro) the PN series appears 
to converge to the "exact" GSF result, (iii) This convergence does not appear to be "monotonic", in the sense that 
the partial PN sums alternate between "overshooting" and "undershooting" the GSF result. [This can be contrasted 
with the case of u l (x), where a monotonic convergence is observed at least through 7PN [3,[l^].] (iv) Close to the 
ISCO the PN series does not show a good convergence; going to higher-order PN does not necessary improve the 
accuracy of the PN approximation there, (v) However, as expected, the PN series does extremely well at small x; for 
x < 1/50 the 5PN model approximates p to within a few parts in 10 4 , even when the (yet unknown) 4PN and 5PN 
constant terms are neglected. 

In the next section we will take a closer look at the GSF data, and explore more quantitatively what can be learned 
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ro/M 


x = M /r 


P 


80 


0.0125 


0.0024117(9) 


57.142. . . 


0.0175 


0.0048913(6) 


50 


0.0200 


0.006494(2) 


44.444. . . 


0.0225 


0.008351(2) 


40 


0.0250 


0.010470(1) 


36.363. . . 


0.0275 


0.0128610(8) 


34.2857 


0.0291. . . 


0.0146099(8) 


30 


0.0333. . . 


0.0195438(4) 


25 


0.0400 


0.0291863(3) 


20 


0.0500 


0.0479916(5) 


19 


0.0526. . . 


0.053862(3) 


18 


0.0555. . . 


0.060857(2) 


17 


0.0588. . . 


0.069279(4) 


16 


0.0625 


0.079537(2) 


15 


0.0666. . . 


0.092199(3) 


14 


0.0714. . . 


0.108061(2) 


13 


0.0769. . . 


0.128280(3) 


12 


0.0833. . . 


0.154578(3) 


11 


0.0909. . . 


0.189605(3) 


10 


0.100 


0.237610(4) 


9 


0.111. . . 


0.305750(5) 


8.5 


0.117... 


0.351000(6) 


8 


0.125 


0.406767(6) 


7.5 


0.133... 


0.47651(1) 


7.4 


0.135... 


0.492527(7) 


7 


0.142. . . 


0.56528(1) 


6.8 


0.147... 


0.607693(9) 


6.5 


0.153... 


0.68059(1) 


6 


0.166... 


0.83413(4) 



TABLE II: Numerical GSF data for p(x). Figures in brackets are rough estimates of the absolute error in the last displayed 
digit (so, for example, the value in the first line stands for 0.0024117 ± 0.0000009). These estimates include the error from the 
finite numerical mesh size (including error from the large I tail approximation [§[), added in quadrature to an estimate of the 
error from residual spurious initial waves. The value for ro = 6M is obtained through an extrapolation based on 18 data points 
(not shown here) between ro = 6.0005M and ro = 6.05M. 

from comparing it with the EOB predictions. 

IV. COMPARISON OF GSF RESULTS WITH EOB/PN PREDICTIONS 

A. Preliminary tests 

We begin by checking that the numerical GSF results are consistent with the EOB/PN prediction (f3"4"|) . with the 
specific analytic values of the coefficients given in Eqs. (|4Tl) and (|42[) . To this end, let us consider the 4PN residual 
quantities 

A 4 (x) ee (p(x) - p 2 x 2 - p 3 x 3 ) /x 4 , 

A%(x) ee (p(x) - p 2 x 2 - p 3 x 3 - p l £ s x 4 bur J /x 4 , 

A++(a;) ee (p(x) - P2 x 2 - p 3 x 3 - p l ° 6 x 4 Inx- p l ° e x 5 ]nx) /x 4 , (44) 
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(log only) 



2PN 



0.02 0.04 0.06 0.08 

x 



0.1 



0.12 0.14 



0.16 



FIG. 1: Numerical GSF data for p{x) compared with various EOB/PN approximations. For clarity we show p{x)/x 2 rather 
than p(x) itself, recalling the small- x asymptotic behavior p(x) oc x 2 . Dots represent the GSF data points shown in Table HT1 
and the blue (darkest) line is a simple cubic spline interpolation of these data. Numerical error bars are too small to show on 
this scale. Other lines display various analytic PN models p PN constructed from Eq. ((34}; the plots labelled 'nPN' show p PN 
through 0(x n ), where at 4PN and 5PN only the known, logarithmic terms are included. The inset shows an expansion of the 
small- a; portion of the plot. Recall x is the (dimensionless) gravitational potential, so x — > corresponds to rn — s* oo. The 
a;-domain here extends to the ISCO location at x = 1/6. 




FIG. 2: The fractional difference jp PN — p\ / p between the various PN models shown in Figure[T]and the numerical GSF data, 
as a function of x. (As in Figure [T] the 4PN and 5PN models include the logarithmic terms only.) Note the logarithmic scale 
of the vertical axis. The magnitude of the "kink" at the far left end of the 5PN plot is well within the numerical error for the 
left-most data point (at x = 1/80). 



where the coefficients p2, P3, p£ s and p£ s are those given in Eqs. (|4Tj) and (|42)l above. Using the notation introduced 
in Eq. (1341) above, we can write the small-x expansion of these quantities as 

A 4 (x) = pl + p 1 ° e \nx + {pt+p l ° s \nx)x + O{x 2+0 ), (45) 
Aj(z) = p\ + (Pi + p l s s \nx)x + 0(x 2+ °), (46) 
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A++(x) = pl + pix + 0(x 2 +°). (47) 

Note that the residue A4 diverges logarithmically at x — > 0, while the "improved" residue A|" is finite at this limit but 
has a divergent derivative there. The "further improved" residue A^~ + , however, is both continuous and differentiablc 
at x — + , with A^~ + (0) = p\ and (A^ + )'(0) = p\. In particular, the task of extracting the unknown PN coefficients 
p\ and p\ from the numerical data amounts to resolving the values of A^~ + (a;) and its derivative at x = 0. 

In Figure [3] we have plotted the functions A^a:), A^(a:) and A^~ + (x) based on the GSF numerical data. A visual 
inspection of the plot reveals the following, (i) The data for A^a:) is consistent with the expected logarithmic 
divergence at x — > 0. (ii) This divergence seems to be eliminated in A^(x), again as expected. A closer examination 
(see the inset) suggests that the derivative (Aj )' begins to vary rapidly at x < 0.05; although we cannot be conclusive, 
this behavior is what one would expect if A4" exhibited a oc x In x term as in Eq. (l46l) . (iii) The data for the residue 
A^ + seems to vary more smoothly as a function of x, even at small x; extrapolating "by eye" to x — > suggests that 
both A^~ + and its derivative attain finite values at x = 0. 

These results give a firmer basis to the conclusion that the GSF data is indeed consistent with the PN expression 
with the specific analytical values of the 2PN and 3PN coefficients given in Eqs. (|41[) . In particular, our results 
confirm the occurrence of a high level of cancellation in the small- a; expansion of the GSF expression for p{x), Eq. 
(1431) . resulting in a oc a; 2 leading-order behavior in agreement with the prediction of PN theory [42|. Furthermore, 
Figure [3] suggests that the GSF data correctly capture the logarithmic terms occurring at 4PN and 5PN, which were 
recently determined analytically. 



200 




x 



FIG. 3: Numerical GSF data for the 4PN residual quantities /S.a{x), A^{x) and A^ + (x) [see Eqs. (14"4")l for definitions]. Curves 
are simple cubic-spline interpolations of the numerical data, and the inset displays the A^ data at a modified aspect ratio for 
better clarity. Error bars are computed from the estimated numerical errors in p(x) indicated in Table [TTJ As discussed in 
the text, these graphs illustrate the consistency of the numerical GSF data with the EOB/PN prediction through 3PN order. 
Furthermore, they suggest that the GSF shows the correct (analytically predicted) logarithmic running at 4PN and 5PN. 

Figure [3] also visually illustrates the difficulties inherent in our comparison. Even though the fractional numerical 
error in p{x) has been kept roughly uniform (^-independent) in our computation, the corresponding error bars on 
the various A4 residues increase significantly with decreasing x (as expected). This, of course, reduces our ability 
to make precise statements about the behavior of the GSF data at small x, particularly at sub-leading PN orders. 
The problem is exacerbated by the aforementioned numerical cancellation of the O(x ) and O(x) terms in Eq. (|43|) . 
which effectively acts to amplify the numerical error in p(x) at small x. To offset this problem it would be desirable, 
in principle, to adjust the accuracy goal set for p(x) as a function of x, so that higher accuracy is sought for smaller 
x. However, this cannot be achieved easily with the current GSF code architecture: Recall that our computational 
cost increases quickly with decreasing a; as a result of the longer dissipation time of the initial spurious radiation. 
Likewise, it would be desirable to obtain accurate GSF data for x values smaller than 1/80, which would give us a 
better handle on the asymptotic PN behavior. This too, however, is extremely difficult to achieve with our current 
code for the same reason as above. 



13 



fit model 


fixed params. 


P2 (best fit) 


X 2 / DoF 


L°°-norm 


p 2PN 


none 




21.5941 


6.8 x 10 7 


2.3 x 10" 1 


p 3PN 


none 




14.5748 


5810 


3.6 x 10" 3 


p 4PN 


none 




14.5135 


5264 


4.8 x 10" 3 


p 4PN+ 


none 




13.9665 


29.4 


6.0 x 10" 4 


p 5PN 


none 




14.0544 


4.08 


2.0 x 10 -4 


p 5PN+ 


none 




13.9721 


0.74 


4.5 x 10" 5 


p 6PN 


none 




14.0106 


0.59 


1.6 x 10" 5 


p 6PN+ 


none 




13.9619 


0.58 


1.7 x 10" 5 


p 7PN 


none 




13.9527 


0.61 


1.7 x 10" 5 


p 7PN 
p 7PN 
p 7PN 


P3 

log; 
P3, Pi 

log; 

P3, Pi, 


P5° g 


13.9946 
14.0015 
14.00002 


0.58 
0.56 
0.55 


1.7 x 10" 5 
1.6 x 10 -5 
1.6 x 10" 5 



TABLE III: Numerical determination of the 2PN coefficient pi from the GSF p(x) data. Each line describes best-fit results for a 
particular PN fit model p nPN or p nPN + [see Eq. (|48[) for notation]. In the upper part of the table we have left all PN parameters 
of the various models as freely-specifiable fit parameters, whereas in the last 3 lines some of the parameters (specified in the 
second column) have been fixed at their known analytic values. The third column of the table displays the best-fit value of pi 
for each model (recall that the analytic value of pi is precisely 14), and the fourth column shows the corresponding value of 
the minimized \ 2 divided by the number of Degrees of Freedom. In the last column we indicate the magnitude of the largest 
absolute difference Ap| between a data point and the value predicted for that point from the best-fit model. 



B. Quantitative test of the GSF data against known EOB/PN terms 

Next, let us inspect the agreement between the GSF results and the analytic predictions of EOB/PN in a more 
quantitative fashion. We will do so by fitting PN models to the numerical p(x) data and assessing the goodness of 
the fits using a standard \ 2 analysis. To avoid the loss of accuracy inherent in multi-dimensional fits [431 ]. we will 
be using a "marginalization" procedure whereby we fit only one PN parameter at a time, then fix the value of that 
parameter (in accordance with the analytic PN prediction) and proceed to inspect the next PN level. At each stage 
of this procedure we will fit our data against a few different PN models. For easy reference, we introduce the notation 

p nPN . ^pn U p £ Q 0(3;™) inclusive, excluding the nPN log term, 
p nPN + : p PN up to 0(x™+°) inclusive, i.e., including the nPN log term, (48) 

where "nPN log term" refers to a term of the form p„ s x n lnx. Our convention is that both p nPN and p nPN + include 
all logarithmic terms occurring at orders 4PN through to (n — 1)PN. 

We begin with the leading 2PN term, p ~ P2X 2 , to check to what extent our GSF numerical data confirm the 
analytically predicted value of the 2PN coefficient: p^ na[ytlc = 14 \y e first consider a set of PN models p 2PN , . . . , p 7PN 
and p 4PN+ , . . . ,/9 7PN+ with all PN coefficients taken as free fit parameters. We least-squares-fit each model to our 
numerical GSF data, weighed by the estimated numerical error from Table|TI](we do that in practice by minimizing the 
corresponding \ 2 function over all fit parameters, using the built-in Minimize[] function available in Mathematica) . 
For each fit model we record the best- fit value of pi, ignoring all other fit parameters for now. We also record the 
value of (the minimized) \ 2 , the value of the L°°-norm (i.e., the maximum magnitude of the absolute difference \Ap\ 
between a data point and the best- fit model), and the number of degrees of freedom for the fit (DoF; the number of 
numerical data points less the number of fit parameters). The results are presented in the upper part of Table Hill 
At the bottom (last 3 lines) of Table Hill we present similar fit results for model p 7PN , where this time we fix some of 
the higher-order parameters P3, P4° s , and p^ s at their analytically-known values. 

There are two important caveats one must bare in mind when analyzing the data in Table [TTT1 First, even if the GSF 
data were exact (i.e., containing no numerical error), still a fit to a high-order PN-like model would not be expected to 
produce the exact analytic values of the PN parameters (p2 — 14, etc.) arising from the systematic, mathematically 
well-defined PN expansion, because the PN series represents an asymptotic expansion while the GSF data includes 
strong-field information. (A corollary is that the "best fit" PN-like model is not at all guaranteed to coincide with 
the actual, mathematical PN expansion.) The second caveat is that, even if the PN model gave a precise description 
of p for all x, still the sum of squares of weighted differences between the GSF data and the PN model would not 
necessarily follow a \ 2 distribution with mean=DoF, since our numerical data are unlikely to represent independent 
sample points drawn randomly from a Gaussian distribution. [Our estimated numerical errors, which we treat as 
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fit model 


fixed params. 


P3 (best fit) 


X 2 / DoF 


L°°-norm 


p 3PN 


P2 


97.953 


3.7 x 10 5 


8.2 x 10" 3 


p 4PN 


(>2 


106.936 


4.9 x 10 4 


1.2 x 10" 2 


p 4PN+ 


t>2 


122.458 


20.5 


4.4 x 10" 4 


p 5PN 


Pi 


120.962 


12.0 


3.6 x 10~ 4 


p 5PN+ 


P2 


124.365 


1.04 


7.8 x 10~ 5 


p 6PN 


P2 


122.256 


0.57 


1.6 x 10" 5 


p 6PN+ 


P2 


123.758 


0.57 


1.6 x 10~ 5 


p 7PN 


P2 


120.914 


0.58 


1.7 x 10" 5 


p 7PN 


loe 

P2, P4 8 


122.929 


0.56 


1.7 x 10" 5 


p 7PN 


log log 
P2, Pi, P 5 B 


122.623 


0.55 


1.6 x 10" 5 



TABLE IV: Numerical determination of the 3PN coefficient pz from the GSF p(x) data. The structure of the table is similar 
to that of Table ITTT1 but here we present fit results for pz (recall the analytic value of pz is 122.6274 . . .). In all the models 
presented we have fixed p2 = 14. The last two lines show results from fits where, in addition, we have fixed the higher-order 
parameter p 4 ° s , and then both p 4 ° s and p]? s , at their known analytic values. 



random statistical errors for the present exercise, may well have a dominant systematic component (particularly from 
residual non-stationarity); furthermore, the numerical errors for different x may well be correlated.] Nonetheless, for 
lack of better options we invoke here the standard x 2 test as a rough measure of goodness for our fits, and accompany 
this with L°°-norm information for a fuller picture. 

With the above caveats in mind, let us inspect the data in Table IIIII We make the following observations, (i) 
The fitted value of P2 matches the theoretical PN prediction to within < 0.4%, so long as the fit model includes all 
terms through 4PN+ at least, (ii) Anecdotically, this agreement becomes extremely good if we allow ourselves to 
fix some (or all) of the higher-order parameters known analytically, (iii) For models p 5PN and beyond the value of 
X 2 /DoF is of order unity, indicative of a statistically good fit throughout the entire ^-domain, and also suggesting 
that our estimates of the statistical error in p are quite reasonable. The slightly low values of x 2 /DoF~ 0.6 (instead 
X 2 /DoF^ 1) probably reflect our somewhat conservative approach to error estimation in our code (see Ref. Q for 
details; for instance, we combine the numerical errors from the various multipole modes by simply adding up their 
absolute magnitudes rather than adding them in quadrature), (iv) The agreement between the fitted pi and its 
analytic value does not seem to improve (or worsen) significantly upon adding fitting terms beyond 4PN+; likewise, 
the value of both x 2 /DoF and the L°°-norm seem to "saturate" beyond 5PN. It appears that higher-order terms have 
magnitudes which are too small (relative to the numerical error) to affect our x 2 analysis (cf. the next subsection and 
Fig. [U where we analyze the magnitude of the various PN "signals" as compared with the "noise" from numerical 
error), (v) Strikingly, the L°°-norm associated with some of the models in Table Hill is extremely low. If one sought 
to model the GSF-computed p(x) at accuracy of (say) better than 10 -4 , then our 5PN+ model would be perfectly 
adequate, even in the strong-field regime. Note that L°° ~ 10 -5 is comparable with the magnitude of the statistical 
error in the GSF data, which explains the saturation of the L°°-norm at this level. 

Having established that the GSF data are quantitatively consistent with the analytic 2PN result p^ nal y tlc = \^ we 
henceforth fix P2 = 14, and turn to consider to what extent the GSF data confirm the analytically predicted value of 
the 3PN coefficient: ^ nal y tic = 122.6274 . . .. Once again, we fit various PN models to the numerical data (keeping 
P2 fixed) and this time record the best- fit values obtained for p 3 . The results are presented in Table [IVJ As before, 
models ( o 4PN + and beyond yield best-fit p% values in a good agreement with the analytic prediction. The level of 
agreement (< 2% difference) is slightly worse than at 2PN; this is expected given the weaker amplitude of the 3PN 
"signal" (cf. Fig. [3] below). Once again, models p 5PN + and beyond yield x 2 /Dof values of order unity and a L°°-norm 
smaller than 10 -4 . Here, too, fixing the known higher-order coefficients improves the agreement with the analytic 
prediction significantly. 

We next apply the same procedure for each of the logarithmic coefficients /9 4 ° s and p l £ s in turn: For the former 
we set both coefficients P2 and ps to their exact analytic values, and for the latter we also fix p 4 ° s . The results are 
presented in Tables |V] and IVI1 respectively. In the case of p 4 ° s the agreement with the analytic value is much less 
impressive than at 2PN and 3PN. Even though models p 5PN+ and beyond again show x 2 /Dof values of order unity, 
the resulting best-fit values for p l £ s vary within as much as ~ 35% of the analytic prediction, and the inclusion of 
additional PN terms does not seem to lead to any convergence of these values. (Note, however, the rather excellent 
agreement when fixing the higher-order logarithmic coefficient P5° s .) The agreement is even less striking in the case 
of pi; 08 (Table |VI|, where the best- fit values for models with x 2 /DoF of order unity vary within ~ 50% of the analytic 
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fit model 


fixed params. 


p 4 og (best fit) 


X 2 / DoF 


L°°-norm 


p4PN+ 


P2, Pi 


177.667 


50.04 


6.9 x 10~ 4 


p 5PN 


P2, Pi 


178.371 


51.39 


7.4 x 10 -4 


p 5PN+ 


P2, Pi 


226.247 


1.54 


3.2 x 10" 4 


p 6PN 


P2, Pi 


169.686 


0.56 


2.0 x 10~ 5 


p6PN+ 


P2, Pi 


194.132 


0.56 


1.6 x 10~ 5 


p 7PN 


P2, Pi 


108.606 


0.56 


1.7 x 10" 5 


p 7PN 


log 

P2; Pi, P 5 


168.474 


0.55 


1.6 x 10" 5 



TABLE V: Numerical determination of the 4PN logarithmic coefficient p 4 og from the GSF p(x) data. The structure of the table 
is similar to that of Tables fTTTl and [TVl but here we present fit results for p 4 og (recall the analytic value of p 4 og is 167.4666 . . .). 
In all cases we have fixed p2 and pi at their known analytic values. The last line shows results from a fit where, in addition, 
we have fixed the higher-order parameter p!° g . 



fit model 


fixed params. 


p' 5 og (best fit) 


X 2 /DoF 


L°°-norm 


p 5PN + 
p 6PN 

p6PN+ 
p 7PN 


log 

P2, Pi, P 4 

loe 

P2, Pi, P 4 

P2, Pi, P 4 ° S 
log; 

P2, Pi, Pi 


-106.14 
-823.96 
-849.15 
-1738.65 


74.33 
0.54 
0.57 
0.55 


7.9 x 10" 4 
1.6 x 10 -5 
1.8 x 10~ 5 
1.6 x 10" 5 



TABLE VI: Numerical determination of the 5PN logarithmic coefficient p 5 og from the GSF p(x) data. The structure of the 
table is similar to that of Tables Mil IIVI and [V] but here we present fit results for p!° g (recall the analytic value of pj° g is 
— 1619.4 . . .). In all cases we have fixed p2, Pi and p 4 ° g at their known analytic values. 

prediction (but note the rather good — perhaps "incidental" — agreement for p 7PN ). 

It is of no surprise that the higher-order PN coefficients are less well determined than the lower-order ones. We expect 
higher-order PN terms to have lower relative magnitudes (or "signal"), especially at small x, and we heuristically 
expect the statistical errors in determining each of these terms to be roughly inversely proportional to its typical 
"signal-to- noise" ratio (with the "noise" here being provided by the statistical numerical error in p). In addition, 
in the 4PN and 5PN cases, the signals from the logarithmic terms are strongly correlated with the signals from the 
unknown 4PN and 5PN 'constant' terms, which further reduces the extraction accuracy of the logarithmic coefficients. 

It is instructive to examine the magnitudes of the various PN signals relative to the numerical noise (as was also 
done in Ref. [15]), and we do so in the next subsection. 

C. "Signal-to- noise" analysis 

We are interested here in a rough estimation of the magnitudes of the various PN contributions to p(x), in comparison 
with the numerical error in p. This requires knowledge of some of the yet-unknown high-order PN coefficients. In 
the next section we attempt to determine some of these coefficients through a systematic analysis. For our present 
purpose, however, we shall content ourselves with a cruder approach, whereby we fit all unknown PN coefficients 
simultaneously using a simple, tentative high-PN model. Specifically, we fit the numerical data to an 8PN model 
in which all known parameters (p2, P3, p 4 ° s and p^ s ) are pre-fixed at their analytic values, and where p\. . . p%, as 
well as pg° s , are left as fitting parameters. Since the form of the logarithmic dependence at high PN order is not yet 
clear [there may well occur cx (lnx) 2 terms, for example] we choose to crudely absorb all possible 7PN terms in a 
term of the standard form p^x , and similarly for 8PN. A least-squares fit then yields the following values (which one 
should consider merely indicative): p\ = 68.48, p% = -4742.81, pg = -771.41, p l ° s = -7349.30, p c 7 = 5757.52, and 
p§ = 11278.11. We use these values to deduce the amplitudes of the PN contributions through 8PN, as functions of 
x. We show these amplitudes in Figure 21 along with the amplitude of the "noise" from numerical error. 

The accuracy Ap n with which we can extract the value of a given PN coefficient p n (assuming all lower-order 
coefficients are known) is roughly inversely proportional to the typical signal-to-noise ratio (SNR) associated with the 
corresponding (individual) nPN signal. Inspection of Figure 0] reveals that the SNRs of the 2PN and 3PN signals are 
roughly equal in the strong field (x > 0.1), but in the weak field the 2PN SNR increases gradually with respect to the 
3PN one, up to about a relative factor 10 at x = 1/80. This is consistent with our finding that Ap 2 is about 5 times 
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FIG. 4: "Signal" amplitudes from various PN contributions to p(x) (thin lines), compared with the "noise" amplitude from 
numerical error (thick red line). For this plot, all unknown PN coefficients through 8PN were crudely estimated by fitting the 
numerical GSF data to a tentative 8PN model. Each of the signal lines displays the amplitude of the total contribution from a 
particular PN order, which at 4PN through 6PN also includes a logarithmic term of the form p n x n lnx (at 7PN and 8PN we 
crudely absorb all possible logarithmic running in an effective term of the form p n x n ). Note the logarithmic scale of the vertical 
axis. The total 5PN contribution changes its sign around x = 0.0535, where the (negative) p%x 5 term conspires to cancel out 
the (positive) p l ^x 5 \nx term. (Note that in this figure the labels 'nPN' indicate the contribution of an individual PN order, 
unlike elsewhere in the text where nPN stands for a partial PN sum.) 



smaller than Ap^. 

The 3PN SNR, in turn, is always larger than the 4PN SNR, by an amount varying from a factor ~ 3 near the ISCO 
to a factor ~ 10 in the weak field. The extraction error Ap l £ s was found to be ~ 15-20 times larger than A/93, which 
is a bit more than one might expect based on the simple T/SNR' scaling argument. However, here one should recall 
that the problem of fitting p£ g also involved fitting the unknown parameter p|, which necessarily has the effect of 
further increasing the error Apjj 08 . Hence, our "empirical" level of uncertainty in p l £ s is not unreasonable. 

The case of the 5PN signal is more involved. As manifest in Figure |H the 5PN signal is strongly suppressed, and it 
entirely vanishes around x = 0.0535. This is a result of a cancellation between the negative term p|x 5 = — 4742. 8I2 5 
and the positive term p£ s x 5 lnx = — 1619.428571x 5 lnx (recall lna; < in the relevant domain). Inspection of the 
separate SNRs from the two 5PN contributions (not shown in the plot) reveals that they are each about equal to the 
4PN SNR in the strong field (x > 0.1), and gradually decrease in the weak field to about 1/10 of the 4PN SNR at 
x = 1/80. However, our simple '1/SNR' scaling argument clearly does not apply here, since the two 5PN signals are 
manifestly strongly correlated between themselves, and, moreover, the effective noise level for the 5PN signals has 
a large contribution from the poorly resolved coefficient p\. These complications suggest that p ° s should be rather 
poorly resolved using our current GSF data, as we indeed concluded above through experiment. [44| 

The information in Figure [4] serves not only to check that our fitting accuracies for the known PN parameter make 
sense, but also (more importantly) to establish an expectation as to which of the yet-unknown PN coefficients we 
might hope to resolve using our current GSF data. It is immediately clear, for example, that we are not likely to be 
able to extract the 7PN and 8PN coefficients at any level of confidence, since the crucial weak-field portion of their 
signal is deeply buried in the numerical noise (and, furthermore, the strong-field portions of the 7PN and 8PN signals 
are likely to be strongly correlated, as is visually evident from the near-proportionality of the two signals in Figure 
E}. 

The situation with the 6PN signal is somewhat better: We have statistically-significant data available for x > 0.02, 
so one might hope to be able to extract some information about the 6PN coefficients. However, since the 4PN and 
5PN parameters are not known analytically, the extraction errors from these parameters will effectively add to the 
noise level relevant for the 6PN signal. Furthermore, since at present the 6PN logarithmic dependence is not known 
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analytically, it might prove difficult to disentangle (with statistical confidence) the pgX e term from any possible 6PN 
logarithmic terms. 

As for the 4PN and 5PN terms: Here the signals lie above the numerical noise across the entire x domain of 
the numerical data, and, in addition, the logarithmic running terms are known analytically. This suggests that our 
numerical data is sufficiently accurate to allow fitting the values of p| and p§. Here too, however, we note that the 
effective noise level for the 5PN term will be increased by the statistical extraction error from the 4PN term. 

The above is, of course, only a heuristic discussion based on the tentative signal amplitudes shown in Figure UJ In 
the next section we carry out a more systematic analysis, aimed at extracting the values of (and determining the error 
bars on) as many PN parameters as is possible given the numerical data. 

V. DETERMINATION OF UNKNOWN EOB/PN PARAMETERS 

The main potential payoff from a GSF-EOB synergy lies in the determination of the strong-field behavior of some 
0(q) functions relevant to the dynamics of binary systems. In this section, however, we remain for the time being 
within the context of PN theory, and attempt to use the GSF data to determine some of the unknown higher-order 
weak-field parameters entering the PN expansion of p(x). In the following section we will then turn to explore the 
strong-field information contained in our GSF data. 

A. Numerical determination of yet-unknown, higher-order PN expansion parameters of p(x) 

The heuristic discussion in the previous section suggested that we might be able to determine p\ and possibly p|, 
while the ability to determine the 6PN parameters remains less certain. The tentative signal amplitudes in Figure 0] 
also suggest that in fitting a PN model to the GSF data we should include all terms down to 7PN, or perhaps 8PN 
(otherwise, the remaining "unmodelled" piece in the numerical data would significantly affect the estimates of the 
lower PN parameters). As a final preliminary note, we observe (extrapolating the upper two curves in Figure [3] "by 
eye" to x = 0) that the value of p\ is expected to fall in the range 50 < p\ < 100. 

Our analysis proceeds as follows. Fixing all known parameters (p2, P3, p 1 ^ 5 , p 5 ° g ) at their analytic values, we 
consider a set of PN models from p 5PN + through to p 8PN . At 7PN and 8PN we attempt a variety of 7PN logarithmic 
dependences, including x 7 lnx and a; 7 (lnx) 2 (each in separate or combined forms). We fit each model to the GSF 
data and record the best-fit values of p|, p§, p% and pg° s , along with the values of % 2 /DoF and the L°°-norm. The 
results are displayed in Table IVIII 

We observe that, in accordance with our expectation, the goodness of the fit (as measured by % 2 /DoF and the 
L°°-norm) "saturates" at the 7PN level; lower-order PN models do not fit the numerical data as well, and the fit does 
not seem to improve when including 8PN terms. We therefore focus our attention on the 7PN and 8PN models in 
Table IVHl These models predict p\ values between 65.6983 and 75.6533. We consider this interval a measure of the 
uncertainty, Ap|, in our determination of p\. We note, however, that many of the p\ values in the table are clustered 
around 68 or 69. We shall take p\ = 69 as our "best guess" value, and as a rough error margin take the asymmetric 
range between 65 and 76. The values predicted for p\ vary between —4444.37 and —5933.41, with many of the values 
clustered around —4700 or —4800. We shall estimate p\ = —4800, with asymmetric error range between —6000 and 
—4400. The values for p£ are clearly dominated by random noise, and provide us with no meaningful information. 
The best fit values p 6 ° s are likewise very "noisy" , although is seems safe to conclude that the actual value of p 6 og is 
negative (and perhaps of order a few thousands). As expected, the values of the various 7PN and 8PN coefficients 
(not shown in Table PVTl} are entirely dominated by noise and cannot be determined. 

In summary, the information included in the currently available GSF data allows us to conclude 

P2 = 69±L Pl = -4800±?™o> P6° S <0. (49) 

We note that the rather large uncertainties in p\ and p\ are in fact similar in magnitude to the ones obtained in the 
previous section for (correspondingly) p^ 6 and p 5 og — cf. Tables IVl and I VII Somewhat disappointingly, the accuracy of 
our numerical data does not allow us to set tighter constraints on pi and pg , nor does it give us reliable access to the 
values of p% and p 6 og . 
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replacing the term p 7 x 7 with p l ° E x 7 In x 
replacing the term p^x 7 with p° s x 7 (\nx) 2 
adding a term p° s x 7 (lnx) 2 
omitting the term p l ^x 7 In x 



TABLE VII: Numerical determination of yet-unknown, higher-order PN expansion parameters of p(x) from GSF data. Each 
line in the table corresponds to a particular PN fit model as indicated in the first column [referring to Eq. (|48[) for notation]. 
In all models we have fixed p2, P3, p^ s and pj° s in accordance with their known analytic values [Eqs. (|41[) and (I42p ]. For each 
model we give the best-bit values of p%, pi, pi and pj° g , and the corresponding values of x 2 /DoF and the L°°-norm. Note the 
last model in the table is the tentative one considered in Subsec. [IV CI 



B. Implications of p\ and p\ for EOB theory 

What can we learn from the above estimates of p\ and p\ about the values of the yet-undetermined 0(q) (logarithmi- 
cally running) EOB parameters a^(\nx), ag(lnx), . . . and ^(lnx), ds(lna;), . . .? To discuss this issue we need to come 
back to the exact EOB relation between the GSF-determinable function p(x) and the basic functions a(x) and d(x) 
of the EOB formalism, i.e., to Eqs. (|27|) -(p0 | . Above, we considered the first two terms in the PN expansion of these 
equations, corresponding to the 2PN and 3PN levels. If we now consider the higher-order terms in the PN expansion 
of these equations, and if we separate them according to whether they contain lnx or not, we find that p\ and p§ are 
related to the coefficients entering the PN expansions of the a{x) and d{x) functions [setting as(lnx) = a§ + a l £ s lnx, 
etc.] via 

27 _ _ Q , 

p\ = — r -7a 4 + 10a c 5 -6d 3 + d c 4 + -a l ° 6 , (50) 



4 u * 2 

— - 14a c 5 + 15a c 6 - U% + d% - 8a!° s + - 



Pi = -— -14a c 5 + 15o§-6d£ + d§ -84 os + -<4 og . (51) 



Substituting from Eqs. ([M]). ([551) and (02) w e then obtain the constraints 

10a^ + ^ + ^4 os ~ 518-6±2, (52) 
14^ + 6^ - 15^ " ^5 + 84° S - y4° S ^ 4779;f™. (53) 

We note the fortunate fact that the large relative uncertainty in our fit for p\ manifests itself only with a modest 
fractional error on the right-hand side of the constraint equation (|5"2"|) [this is because the term p\ ~ 69 in Eq. (|5"0|) 
happens to be quite a bit smaller than the known part 7a 4 + 6d 3 ~ 442.815 in that equation]. On the other hand, the 
fractional uncertainty on the right-hand side of the constraint equation (|53p remains rather large. 

As noted in Ref . [28| , knowledge of the PN expansion coefficients of p(x) does not on its own allow us to constrain 
separately the PN expansion coefficients of the EOB functions a(x) and d(x), but only combinations thereof. Another 
subtlety which enters our constraints (|5"2"]) , (f53"|) is that their left-hand sides involve a combination of non- logarithmic 
and logarithmic PN coefficients. Our GSF results do not either give us a direct access to the separate values of a£ s 
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and ag° s , but only to the specific combinations of a£ s , d^ 6 , Og° s and o^° s that enter p l £ g and p^ 6 , namely 

P Y = +104 OS + 4 OS , (54) 
= -144 og + 15a 6 os -64 og + 4 os . (55) 

Let us note in passing that Ref. [28| has discussed an approximate way of estimating some 'effective' values of 
the coefficients as and oq, roughly corresponding to averages of the logarithmically running parameters ag(lnu) and 
ae(lnu) over an interval of u around the ISCO. With the present limited information that we can derive from our 
GSF data, we cannot meaningfully combine these approximate estimates with our results (I5"2"1) . ([53")) . Ref. 28] went 
on to suggest gauge-invariant quantities other than p (namely, the "whirl" frequency and angular momentum in a 
zero-binding zoom- whirl orbit), which give a handle on the strong- field behavior of the a function without involving 
the d function. We shall come back below to possible ways of combining the knowledge of such additional quantities 
with the results that can be derived from our present GSF data. 

VI. ON THE DETERMINATION OF THE GLOBAL STRONG-FIELD BEHAVIOR OF VARIOUS 0(u) 

EOB FUNCTIONS 

So far we have used our p(x) data to (i) test the GSF (and the PN) calculation(s) and (ii) constrain some of the 
unknown higher-order PN parameters entering the weak-field limit of EOB theory. However, as already mentioned 
above, the main potential payoff from a synergy between GSF and EOB formalisms is the determination of the strong- 
field behavior of some 0(v) functions relevant to the dynamics of the binary system. Indeed, the EOB formalism has 
shown its ability at accurately describing all phases of the merger process of comparable-mass binaries, from the early 
inspiral to the final ringdown [23, 24], in terms of a few basic functions, and notably the functions A(u, v) and D(u, v) 
whose 0(y) parts are the functions a(u) d(u). Our numerical results for p(x) therefore give us, through the relations 
(f2"T)) - (|3"0)) . a first-ever access to the strong-field behavior of (a combination of) functions entering the description of 
the binary's conservative dynamics. As such, our results have the potential to inform the development of the EOB 
formalism, which is of great interest and timeliness in the context of the vigorous effort to model gravitational-wave 
sources for existing and planned detector projects. In this context, it would be desirable to replace our sample 
of tabulated values of the function p(x), Table [TTl by some analytical representation which faithfully matches our 
numerical results throughout the entire domain < x < 1/6, and which is likely to remain adequate even for larger 
values of x. 



A. Accuracy threshold on the global representation of p(x) 

Before discussing various ways which might be used for obtaining such global analytical representations of our 
strong-field data, one should start by assessing the accuracy requirements that such global representations must 
satisfy. In the following we suggest a way of quantifying the desired accuracy standard for p(x). 

The currently most accurate version of the EOB formalism [23| has found, especially for the equal-mass case, an 
excellent fit (within the NR error bar) between the EOB waveform and the NR Caltech-Cornell waveform [36| all over 
a certain thin 'banana-like' region in the plane of the two effective parameters a c ffcctlvo ; fl cffcctivo w hi cn are use d to 
parametrize the form of the global, Pade-resummed function Ap(u, v\ a| ffectlve ; a effective^ Comparing the variation of 
Ap along the center of the good- fit region (for v — 1/4) with the variation of an A function of the type considered 
here [i.e., A(u, v) — 1 — 2u + va{u) + 0(i/ 2 ), where a(u) is allowed to change while the 0(v 2 ) terms are kept fixed], 
allows us to relate the variation of Ap(u) to the variation of va(u) (with v = 1/4). Dividing by v = 1/4, one can then 
use this variation of the function a(u) to infer the corresponding variation of the function p(u) using the relations 
(|2T)) - ([5ni) . (In these mathematical expressions one can rename the variable u as x.) Using this procedure, one finds 
that the resulting variation of the function p(x), as one moves along the center of the good-fit region, stays globally 
quite small all over the strong-field interval < x < 1/3 explored by the EOB evolution. In particular, restricting 
to the interval < x < 1/6 currently accessible to GSF methods, one finds that the variation of p(l/6) which is 
compatible with the current NR/EOB agreement is between —0.0028 (for the lefmost part of the banana-like good-fit 
region) and +0.0028 (for its rightmost part). [Here, we use as a reference point for evaluating the variation the value 
of p{x) corresponding to the 'intersection values' ( a | ffective = -22.3, a^ ffcctivc = +252) selected in Ref. [Hj].] When 
extending the strong-field interval as far as the 'light ring' (45[, i.e., x — 1/3, one finds that the 'allowed' variation of 
p(l/3) is somewhat larger, and of order ±0.01 (and even more in the leftmost region). 

Summarizing: The minimal accuracy requirements that one should impose (at present) on a global representation 
of the function p[x) is a maximum abolute difference (i.e., an L°°-norm) smaller or equal to 2.8 x 10 -3 over the 
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interval < x < 1/6, and staying smaller than about 10 2 over the full interval < x < 1/3 where GSF p-data might 
eventually be available. 

B. On the use of actual, PN expansions for representing p(x) globally 

One can think of several options for constructing global analytical representations of the function p(x) from our 
GSF results. 

A first option might be to try to represent p(x) by a sufficiently large number of terms of its (actual) PN expansion 
in powers of x (including all needed logarithms). However, our Figure [4] clearly shows that this is an impossible task. 
Even if one limited one's ambition to representing p(x) by its PN expansion on the interval < x < 1/6, Figure 
|4] shows that one should first determine the PN expansion coefficients of p(x) beyond the 8PN level, and possibly 
beyond the 9PN level. Indeed, if we use the fitted values of p7 and p% quoted above (where, for simplicity, we drop the 
superscript c) as being indicative of their real values (when absorbing the logarithms in the corresponding PN terms) , 
we have p-jx 7 ~ 0.0206(6a;) 7 , and psx & ~ 0.0067(6cc) 8 , which fail to reach the required level 0.0028 for x = 1/6. From 
the ratio of these two terms (D'Alembert criterion) [46[ one expects the 9PN term to be of order 0.0022(6x) 9 , and to 
barely meet our requirement. However, it is currently unthinkable to be able to either GSF-numerically compute all 
the needed PN coefficients with a decent accuracy [47| or to derive them analytically. Moreover, even if we had at 
hand these higher-order PN terms, they would be essentially useless for controlling the value of the function p{x) in 
the strong-field interval 1/6 < x < 1/3 which is crucial for the applicability of the EOB formalism. Indeed, the formal 
computation of the 7PN , 8PN and 9PN contributions we just wrote down in the extended interval 1/6 < x < 1/3 
shows that the PN expansion completely loses its numerical validity when x ~ 1/3 (for instance psx 8 ~ 1.7(3x) 8 is 
formally of order unity when x ~ 1/3). Evidently, all this is an illustration of the fact that the PN series, being a 
weak-field expansion, cannot be expected to cover the very-strong-field regime 1/6 < x < 1/3 corresponding to radii 
ranging between the ISCO and the light-ring. [Recall, however, that the Pade-resummed functions used in the EOB 
formalism seem adequate to accurately describe this regime.] 

C. On the use of effective, PN-like expansions for representing p(x) globally 

If the actual PN expansion of the function p(x) (corresponding to the mathematically exact Taylor-type expansion 
around x = 0) is inadequate for defining a global representation of this function, one might think that some type of 
effective PN-type expansion of p(x) might do a better job. Indeed, we have seen above that by fitting various high- 
order PN-type models to our numerical GSF data we could represent the function p{x) in the interval < x < 1/6 
by polynomials (with extra logarithms) with L°°-norm values as small as ~ 10~ 5 . Two of the remarkable aspects of 
our 'experimental' findings in Sec. IIVI above were: (i) one could obtain L°°-norms of order ~ fewlO -5 (i.e., at the 
level of our statistical numerical error) by using much fewer terms in the PN-like expansion than expected from our 
rough convergence analysis in the previous subsection (e.g., 5PN with logarithms, instead of 9PN from the expected 
convergence rate); and (ii) one could (nearly) obtain the required accuracy level L°°-norm ~ few 10~ 3 by using as low 
a formal PN accuracy as 3PN (when fitting for p2 and p3 without fixing them by our analytical knowledge). However, 
these results apply only to the restricted interval < x < 1/6, and do not give us any control on the type of PN-like 
expansions needed to globally represent p(x), within the accuracy threshold mentioned above, on the doubled interval 
< x < 1/3. From Weierstrass' Approximation Theorem we know that the (arguably) continuous function p(x) [or 
rather p{x)\ can be uniformly approximated, on any compact interval, by some finite polynomial in x to any required 
accuracy, but we do not know anything about the order of this polynomial. In addition, if we lose the connection 
between the coefficients of this polynomial and the Taylor coefficients of p(x) at x = 0, we lose the advantage of having 
an advance analytic knowledge of some of these coefficients. Lastly, we should recall that constructing an accurate 
PN-like model based on a fit to GSF data entails having at hand a dense sample of GSF data points across the entire 
range < x < 1/3. Obtaining such a data set can be extremely computationally expensive. 

We conclude from this discussion that any type of PN expansion (actual or effective) is ill-suited for obtaining a 
global representation of p(x) that meets the rather strict accuracy requirements discussed above in the full desired 
interval < x < 1/3 (of which only the first half has yet been covered by our current GSF results). 

D. Accurate global representations of p(x) using multiple-point Pade approximants 

We wish to introduce here a new strategy for defining sufficiently accurate global representations of the strong-field 
behavior of dynamically useful functions based on the combined use of a minimal amount of strong-field information, 
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together with the currently available weak- field (i.e., PN) information. This strategy has some similarity with the 
strategy used in defining the EOB basic functions, and notably A(u, v), by means of Pade approximants, but it goes 
further in its way of incorporating some strong-field information. Indeed, the Pade approximants used so far in EOB 
theory have always been one-point Pade approximants (i.e., rational functions constrained to reproduce some given 
Taylor series around, say, the point x = 0). By contrast, the general strategy we propose here consists in using 
multiple-point Pade approximants, i.e., rational functions constrained to reproduce several given Taylor series around 
several points, say x\ = 0, x%, . . .. One of the expansion points at which one injects some (PN) information is still 
the weak-field limit x = 0, but the other points x 2 , %3, ■ ■ ■ will be taken in the strong- field domain, and will be used to 
inject a limited, but crucial amount of strong-field information (obtained, say, by a correspondingly limited number 
of GSF computations). Here we shall demonstrate the applicability of our method by considering the function p(x) 
48]. However, the basic idea is clearly very general, and we believe it can provide a powerful tool for advancing AR 
technology. 

Let us start by considering a very simple "two-point" Pade approximant of p(x), based only on readily available 
p-information given at x — and at the ISCO, x — 1/6. Specifically, let us try to base our approximant on only 
4 "easy" pieces of information, namely the values p"(0), p"'(0), p(l/6) and p'(l/6) (where, recall, a prime denotes 
d/dx). The first two values [which themselves assume the knowledge that p(0) — — p'(0)] are simply related to the 
2PN and 3PN parameters p2 and p3, which are given analytically in Eq. (|4Tj) . The latter two values are obtained from 
the strong-field GSF data: p(l/6) is given in the last line of Table [TTl and to obtain the derivative p'(l/6) we used a 
quadratic fit based on 18 near-ISCO data points between x — 1/6.05 and x = 1/6.0005 [the same data points used to 
obtain p(l/6) itself]. We find p'(l/6) ~ 12.66. We then consider a simple 4-parameter Pade model of the form 

p(x) = c x 2 ( \ +C ] X , 2 ) , (56) 
\1 + dix + d 2 x^ J 

where we have factorized the leading-order PN power x 2 . The 4 model parameters {co, c%, d±, d 2 } are uniquely 
determined by our 4 "data points" . Their values are found to be 

co = Pi = 14, ci = 13.3687, d x = 4.60958, d 2 = -9.47696. (57) 

The resulting 2-point-based Pade model for p(x) is plotted in Figure [5] against the actual numerical data. The 
agreement appears strikingly good throughout the entire x range, despite the fact that our model only uses information 
at the two "end" points x = and x = 1/6. The L°°-norm of the full numerical GSF data, with respect to this model, 
is found to be as small as ~ 2.4 x 10~ 3 . Note that this is within the error tolerance stated above, of < 2.8 x 10~ 3 . It is 
remarkable that our extremely simple model (f56f — which is based merely on 4 readily available pieces of PN and GSF 
information — would appear to be perfectly adequate for our current purpose (at least in the interval < x < 1/6). 

In the above analysis we may, of course, replace the derivative p'(l/6) with the value of p at a different strong-field 
point — say x = 1/8 — which can give us a better handle on the "curvature" of p{x). The quartet {p 2 , P3, /o(l/8), p(l/6)} 
again determines the parameters {cq, Ci, d\, d 2 } in the model (l56l uniquely. This time we find L°° ~ 2.4 x 10~ 4 , 
representing a full order-of-magnitude improvement with respect to the 2-point model considered before. Such a 3- 
point Pade model, which still uses a minimal amount of readily-available GSF data, now qualifies for the development 
of EOB models which are ten times as accurate (this might become necessary in the future) . We may obtain similar 
3-point models replacing p(l/8) with other strong-field data points, and in Table |VHT] we list a few such models, 
showing in each case the values of the model parameters as well as the corresponding L°°-norm (over the interval 
< x < 1/6). Such models yield similar values of L°°-norm (with the best performance apparently achieved taking 
the third data point to be in the vicinity of x = 1/9). To improve the model further requires additional parameters 
and, commensurably, additional data points. The last two lines of Table I Villi present examples of 4- and 6-point 
models. 

The above experiments illustrate an important point: It is possible to obtain an accurate global model for p{x) using 
the available PN expression in conjunction with just a handful of near-ISCO GSF data points; a Pade approximant 
based solely on this information can faithfully "bridge" over the entire x domain where stable circular orbits exist. 
This may hold true for other quantities that characterize the 0{v) dynamics, in which case one has a (relatively) 
computationally cheap way of modelling such quantities, which does not require an exhaustive GSF survey of the x 
domain. Of course, in each case one must have at hand some GSF data at intermediate radii in order to be able to 
confidently assess the faithfulness of one's Pade model. The following general scheme suggests itself: (1) Construct 
a Pade model (or a set of models) based on all available PN information and on a handful of readily computable 
strong-field GSF data points; (2) Compute the GSF for a small, complementary sample of intermediate x values and 
use these data points to assess the performance of the Pade model(s). This scheme could save the need to obtain a 
dense sample of GSF data across the entire x domain (as we have done in Sec. Ill of this work). Furthermore, in this 
procedure the "verification" GSF data (which may be most computationally expensive) need not be very accurate as 
in our present work — their accuracy should only match one's error tolerance. 
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FIG. 5: Global fit for p(x) based on the simple 2-point Pade model (|56p [with Eq. I|57|)]. Thick (blue) points are the numerical 
GSF data, while the solid (magenta) line shows the 2-point Pade model, which is based solely on four pieces of information: 
two at x — (the 2PN and 3PN coefficients) and two at the ISCO [p(l/6) and p'(l/6)]. For comparison, 'broken' curves show 
various analytic PN approximations: 2PN (dotted), 3PN (dashed), 3PN including the 4PN logarithmic term (dash-dot, below 
the numerical data points; light green online) and 3PN including both 4PN and 5PN logarithmic terms (dash-dot, top curve; 
dark blue online). Evidently, a mere knowledge of the GSF (and its derivative) at x = 1/6 already improves significantly our 
ability to construct a faithful model of p in the strong-field regime. 
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TABLE VIII: Simple multiple-point Pade models for p(x). Each line of the table corresponds to a particular model, based only 
on p-information at the x-values indicated in the first column. In all models we have used the known values of the 2PN and 
3PN coefficients (the data point u x — 0"). For all data points i / we have used only the value p(x) [e.g., for x — 1/7 we 
have used p(l/7)], except in the 2-point model presented in the first line, where we have used the derivative p'(l/6) in addition 
to p(l/6). The Pade models have the general form p(x) = 14a; 2 (l + cix + C2X 2 + C3X 3 )/(1 + dix + d2X 2 + d3X 3 ), where for the 
4-point model we set C3 = ^3 = and for the 2,3-point models we additionally set C2 = 0. For each model, the values of all 
non-zero parameters c„ and d n (shown in the table) are uniquely determined by the p-information assumed. The last column 
shows the L°°-norm (in the interval < x < 1/6) of the full numerical GSF data with respect to each of the models. 



The fits we just presented for the function p(x) were mainly given to illustrate a general strategy for an efficient 
synergy between GSF, PN and EOB theories. In Ref. [28( one of us proposed several other gauge-invariant quantities 
that may be accessible to both EOB/PN and GSF treatments and that could be used to further constrain the EOB 
functions at O(v). In particular, it was proposed to utilize two gauge-invariant quantities defined for a marginally- 
bound "zoom-whirl" orbit, namely the GSF-corrected values of the azimuthal "whirl" frequency and of the total 
(dimensionless) orbital angular momentum. It was shown that these two quantities uniquely determine the values of 
the crucial 0(v) EOB radial potential a(x), and of its x-derivative, at the strong- field point x = 1/4, i.e., significantly 
beyond the ISCO. This would allow us to extend our strong-field knowledge beyond the interval < x < 1/6 that 
our present study has been limited to. Most importantly, the application of our general strategy to the function a(x) 
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might give us a good handle on the shape of a(x) in the very-strong-field regime. Indeed, by combining the current 
PN knowledge of a(x) near the weak-field point x = (i.e., the information that <z(0) = a'(0) = a"(0) = 0, together 
with the non-zero values of the two numbers 03 and 04) with the hopefully obtainable GSF knowledge of a(l/4) and 
a'(l/4), we might obtain an accurate global representation of the function a(x) in the entire interval < x < 1/4. 
This would be a most valuable information for the EOB formalism. Then, combining this global knowledge of the 
function a(x) with the direct GSF knowledge of the function p(x) obtained here, we could further derive the global 
shape of the function d(x) (at least in the interval < x < 1/6). 

VII. SUMMARY AND OUTLOOK 

We presented here a calculation of the GSF correction to the precession rate of the orbital periastron for a particle 
of mass m in a slightly eccentric orbit around a Schwarzschild black hole of mass M ^ m. Our calculation is exact 
at first order in the mass ratio q = m/M (within a controlled numerical error of 10~ 4 fractionally). Our first main 
result was the computation of a dimensionless measure of the 0{q) contribution to periastron advance, p, as a function 
of the dimensionless azimuthal-frequency parameter x — [Gc~ 3 (M + 77j)f2 v ] 2 / 3 . We computed the function p(x) for a 
dense sample of x- values in the interval < x < 1/6 corresponding to a sequence of stable circular orbits (down to 
the ISCO, located at x = 1/6). 

The bulk of this paper concerned itself with a comparison of the GSF results with PN predictions formulated via 
the EOB formalism, and we explored what can be learned from such a comparison. We demonstrated three main 
goals that can be achieved by combining GSF and EOB/PN results: (i) test the GSF computation and confirm 
the EOB/PN results (and, indirectly, also reaffirm the regularization procedures underpinning both GSF and PN 
theories); (ii) calibrate yet-unknown high-order parameters in the EOB/PN expansion; and (iii) engineer a faithful 
global model for the dynamical quantity in question, which is valid both in the weak field and in the strong field, and 
anywhere in between. 

The new quantitative results of this work are (i) the 0(q) periastron-advance function p(x) (tabulated in Table ITT1 
and plotted in Figure [lj, (ii) an estimate of the 0(q) parts of a few unknown PN parameters [summarized in Eq. 
(|49[) ]. and (iii) an approximate analytic formula for p(x) [Given in Eq. (|56l) with Eq. (|57[) ; alternative, more accurate, 
analytic models for p(x) are described in Table fVIII| . In addition, we have provided numerical confirmations of 
several analytically determined PN/EOB parameters, namely the following PN expansion coefficients of p(x): the 
2PN coefficient P2, the 3PN coefficient p$, and the recently computed logarithmic 4PN and 5PN coefficients p l £ s and 
P5 6 whose knowledge was quite crucial for determining the 4PN and 5PN non-logarithmic contributions. 

Our best estimates for the unknown (non-logarithmic) 4PN and 5PN parameters (let alone the 6PN parameter) 
are rather crude — see the large error bars in Eq. (14^1) . At fault are two main factors. First, the construction of 
p{x) from the various GSF components, prescribed in Eq. (|43[) . involves a delicate cancellation of the two dominant 
terms in the x expansion around x — 0, namely the O(x ) term and the 0(x) term, leaving a final p(x) of 0(x 2 ); we 
found analytically that F£ ilc and F[ , which are both ex x 2 , cancel in the leading O(x ) term, while the cancellations 
occurring in the next (1PN) O(x) term involve Fl cx x 1 / 2 together with the sub-leading, 1PN corrections to F£ irc and 

F[. These cancellations effectively amplify (in relative terms) the numerical error in the small- X GSF data by 2-3 
orders of magnitude. The upshot is a rather poor accuracy in p{x) for small x values, despite the rather high accuracy 
of the GSF data itself. The second limiting factor is the restricted capacity of our GSF code to deal with relatively 
small- a; (i.e., large radii) orbits. In our current time-domain architecture, the computational cost increases rapidly 
with the orbital radius, as a result of the longer evolution time necessary to eliminate the effect of spurious initial 
radiation. We are currently practically unable to compute the GSF for orbits with radii larger than ~ 100M. Future 
advances in GSF computational technology (e.g., better treatment of initial conditions, mcsh-rcfincmcnt techniques, 
or spectral treatments of the perturbation equations — all of which are presently being studied) are certain to reduce 
the computational cost of GSF calculations, and to allow us to obtain more accurate GSF data and for larger radii. 
It would be interesting to revisit our analysis once more accurate GSF results are at hand. 

Most importantly, our determination of the function p(x) on the interval < x < 1/6 has allowed us to compute for 
the first time the strong-field behavior of a combination of the EOB functions a(u) and d(u), thereby giving us some 
information of direct significance for constructing accurate analytical models of the late stages of the dynamics of binary 
systems. We have also indicated at the end of the last section how the GSF computation of the special marginally- 
bound "zoom-whirl" orbit suggested in Ref. [28[ might, thanks to the multiple-point Pade strategy introduced here, be 
used for deriving the global shape of the crucial O(v) radial potential a{x) in the entire domain < x < 1/4. Then, 
combining this global knowledge of the function a(x) with the direct GSF knowledge of the function p{x) obtained 
here, we could further derive the global shape of the function d(x). This would be a most valuable information for 
the EOB formalism, and could significantly help the development of AR models of coalescing binaries. Unfortunately, 
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our current GSF code cannot handle orbits which are not strictly bound, and we are therefore unable to supply the 
necessary GSF data at the present time. We are, however, investigating ways to make this calculation possible, and 
we are hoping to present such analysis in a future paper. 

Finally, we recall that we now have at hand a GSF code for orbits with arbitrary eccentricities [8j . In a forthcoming 
paper [13| . two of us will present a computation of a certain gauge- invariant relation characterizing the conservative 
0(q) dynamics of such orbits. A GSF-EOB analysis may then, for the first time, allow us access to the strong-field 
behavior of the EOB Q function, which enters the description of the "radial" component of the binary motion. It is 
our impression that there are profound prospects for further fruitful synergy between the GSF and EOB frameworks, 
and we expect this activity to gain considerable momentum in the coming years. 
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[37] Refs. HQB1 did determine the strong-field behavior of the gauge-invariant function ifitly), but this function cannot, 

as far as we know, be used to inform the AR description of the dynamics of binary systems. 
[38] In some of the PN/EOB literature the symmetric mass ratio is denoted rj. 

[39] In principle, GSF methods could explore the function p(x) in the full interval < x < 1/3 where circular geodesic orbits 
exist. However, the present structure of our GSF code makes it impossible to explore the perturbations of the unstable 
circular orbits below the ISCO, i.e., in the range 1/6 < x < 1/3. 

[40] This replacement is only a change in the name of a mathematical argument. It differs from the physical statement that 
u = x + 0{y) along the sequence of circular orbits. 

[41] This is a combined result of (1) initial radiation from the particle having to travel longer before being scattered off the 
strong-field potential barrier surrounding the black hole, and (2) the fact that the amplitude of residual late-time decay 
tails increases with decreasing frequency, as low-frequency radiation is less susceptible to black-hole absorption. 

[42] The 2PN nature of p(x) is directly linked to the fact that the effective EOB metric of a binary system starts to differ from 
the Schwarzschild metric only at 2PN [lq ]. 

[43] If we wanted to check simultaneously the numerical values of the 4 analytically computed PN parameters, p2, P3, p 4 ° s , and 
P5° g , we would need to fit for 6 PN parameters at once, including the two extra non-logarithmic terms at 4PN and 5PN. 

[44] One might attempt to make a more quantitative prediction for Ap!° E through analysis of the variance-covariance matrix 
defined on the signal parameter space. Here, however, we content ourselves with a more heuristic discussion. 

[45] In the EOB formalism the functions A(u), D(u), etc. are supposed to have no singularity at least up to values of their 
arguments of order u ~ 1/2 + 0{v) corresponding to the EOB formal analog of the 'horizon'. But a crucial role is played 
by their behavior up to the EOB analog of the light-ring, i.e., for x = 1/3 + 0(y). Note that Eqs. (|27fl - (|30[) suggest that 
the function p(x) will be singular at x = 1/3 [because of the pe(x) contribution]. However, our discussion is meaningful 
for the 'corrected' function p(x) = p(x) — pe{x) = p a (x) + Pd(x), which is indeed predicted by the EOB formalism to be 
regular up to it ~ 1/2 + 0(v). 

[46] The ratio ps/ pi — 1-96 suggests a radius of convergence around x ~ 1/2; this, in turn, suggests that the exact PN series 
might still converge when x ~ 1/3, but at an extremely slow rate ~ C(2x) n . To be precise one should discuss here the 
convergence of the 'corrected' function p(x) = p(x) — pe{x). Indeed, the piece pe{x) will ultimately limit the convergence 
radius of p(x) to 1/3 because of its singular behavior there. However, one can check that, even at the 9PN level, the 
contribution of the PN expansion of pe(x) to that of p(x) is numerically small compared to p7, ps and our estimated pg, 
and would therefore not affect our conclusion. 

[47] We recall that the extraction of accurate PN coefficients requires accurate GSF data at large binary separations (small a;) , 
and that obtaining such data is extremely computationally costly for our eccentric-orbit problem. 

[48] It might even be better to apply it to the function p(x), but we shall not bother to do that here, because we have in mind 
other applications, which we mention below. 



